MODULE MATVAR
!
! Verwendungszweck:
!
!      Bereitstellung der Materialfunktionen
!
!      VORSICHT:  In der Anwendung auf Scheiben bezieht sich Pgas0 nur auf den Zustand
!                 in der Mittelebene der Scheibe (z=0) Die andere Variablen (Tgas, Opazitaeten)
!                 werden hingegen ueber die Dicke der Scheibe als konstant angenommen.
!
!-----------------------------------------------------------------------
!
!   EOSfile    : Filename of the EOS-File
!   OPArosfile : Filename of the Opacity-File -- Rosseland-mean
!   OPAplkfile : Filename of the Opacity-File -- Planck-mean
!
!   Planck     : Flag for using separate Planck-mean opacities - otherwise OPAplk = OPAros is used
!
!   Pgas0      : Gasdruck zum neuen Zeitpunkt
!   Pgas0Z     : Gasdruck zeitlich zentriert
!
!   dPgas0dE   : Ableitung des Gasdruckes nach der inneren Energie
!   dPgas0dD   : Ableitung des Gasdruckes nach der Dichte
!
!   dPgas0ZdE  : Ableitung des zentrierten Gasdruckes nach der inneren Energie
!   dPgas0ZdD  : Ableitung des zentrierten Gasdruckes nach der Dichte
!
!   Nach dem selben Schema Variablen fuer:
!
!   Tgas       : Gastemperatur
!   OPAros     : gas opacity -- Rosseland-mean
!   OPAplk     : gas opacity -- Planck-mean
!
!-----------------------------------------------------------------------

      use config, only : rkind, np

      implicit none

      private
      public :: MATFCT, UPDATE_MATFCT, &
                EOSfile, OPArosfile, OPAplkfile, Planck, &
                Pgas0, Pgas0Z, dPgas0dE, dPgas0dD, dPgas0ZdE, dPgas0ZdD, &
                Tgas,  TgasZ,  dTgasdE,  dTgasdD,  dTgasZdE,  dTgasZdD, &
                OPAros, OPArosZ, dOPArosdE, dOPArosdD, dOPArosZdE, dOPArosZdD, &
                OPAplk, OPAplkZ, dOPAplkdE, dOPAplkdD, dOPAplkZdE, dOPAplkZdD, &
                gamma, xmol


      character(len=62):: EOSfile, OPArosfile, OPAplkfile
      logical          :: Planck

      real(kind=rkind) :: Pgas0(np), Pgas0Z(np), dPgas0dE(np), dPgas0dD(np), dPgas0ZdE(np), dPgas0ZdD(np)
      real(kind=rkind) :: Tgas(np),  TgasZ(np),  dTgasdE(np),  dTgasdD(np),  dTgasZdE(np),  dTgasZdD(np)
      real(kind=rkind) :: OPAros(np), OPArosZ(np), dOPArosdE(np), dOPArosdD(np), dOPArosZdE(np), dOPArosZdD(np)
      real(kind=rkind) :: OPAplk(np), OPAplkZ(np), dOPAplkdE(np), dOPAplkdD(np), dOPAplkZdE(np), dOPAplkZdD(np)

      real(kind=rkind) :: Pgas0A(np), TgasA(np), OPArosA(np), OPAplkA(np)

      real(kind=rkind) :: gamma
      real(kind=rkind) :: xmol

CONTAINS





SUBROUTINE MATFCT
!
! Verwendungszweck:
!     Berechnung der Materialfunktionen und deren Ableitungen - sowohl fuer den
!     neuen als auch den zentrierten Zeitpunkt
!
!     Alle Materialfunktionen werden fuer die Punkte 1 bis np aufgerufen - also auch
!     fuer die an sich bedeutungslosen skalaren werte mit i = npt
!
!     Die Subroutine EOS berechnet die Ableitungen der Materialfunktionen immer gegen das jeweilige Argument.
!     Daher beziehen sich die Ableitungen z.B. auf die Dichte in der Mittelebene (rho0) und
!     nicht auf die primaere Dichtevariable, die integrierte Saeulendichte.
!
!     Da in der aktuellen Version die zeitlich zentrierten Materialfunktionen durch Mittelung ausgerechnet
!     werden und nicht direkt aus den Tabellen mit XZ als Eingangswert haben wir, z.B., TgasZ(rho,rhoA,e,aA)
!     und nicht TgasZ(rhoZ,eZ). Daher beziehen sich die Ableitungen der zeitlich zentrierten
!     Materialfunktionen auch auf X und nicht auf XZ.
!
!      VORSICHT:  Das Format der Realvariablen in den Materialfunktions-Routinen (definiert in mater/com/mater_module.f90)
!                 muss mit dem des SHD-Codes (in code/module.f90 festgelegt) uebereinstimmen !!!!
!
!      VORSICHT:  Die Zustandsvariablen sind einheitenbehaftet - cgs
!
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   dOPAdTgas    : Ableitung der Opazitaet (zum neuen Zeitpunkt) nach dem Gasdruck
!                  wird fuer Rosseland als auch Planck opacitaeten verwendet
!
! Verwendete Unterprogramme:
!
!   EOSINI(EOSfile) : Einlesen der EOS-Tabelle - "EOSfile" ist der Filename inkl. Pfad
!
!   EOS : Berechnung der Zustandsfunktionen
!
!      EOS (ZUFU, MODE, N, RHO, ENE, Z, DZDR, DZDE)
!
!      mit folgender Bedeutung:
!
!      ZUFU:  Auswahl der Zustandsfunktion
!             Gasdruck "PRE", Gastemperatur "TEM", Nabla-adiabtisch "NAB", spez. Waerme cp "CEP"
!      MODE:  "NEW" neue Druck und Energie-werte,  "OLD" die selben wie beim vorherigen aufrufen.
!      N:     Laenge der Zustandsvektoren: RHO, ENE, Z, DZDR, DZDE
!      RHO:   Dichtewerte
!      ENE:   Energiewerte
!      Z:     zu berechnende Zustandsfunktion
!      DZDR:  deren Ableitung nach der Dichte
!      DZDE:  deren Ableitung nach der Energie
!      TXMAT: Zustandsgleichungskennung  CHARACTER(len=70)
!
!
!   OPAINI(OPAfile) : Einlesen der Opazitaets-Tabelle - "OPAfile" ist der Filename inkl. Pfad
!
!
!   OPACITY: Berechnung der Gasopazitaeten
!
!      OPACITY (N,RHO,TEG,OPG,DOPGD,DOPGT)
!
!      mit folgender Bedeutung:
!
!      N:     Laenge der Zustandsvektoren: RHO, TEG, OPG, DOPGD, DOPGT
!      RHO:   Dichtewerte
!      TEG:   Temperaturwerte
!      OPG:   Gasopazitaet
!      DOPGD: deren Ableitung nach der Dichte
!      DOPGT: deren Ableitung nach der Temperatur - koennen wir so nicht brauchen -> umrechnen in Dichte & Energie
!
!-----------------------------------------------------------------------

      use config,  only : rkind, np
      use primvar, only : X, ME
      use physco,  only : z1
      use zvar,    only : rho0
      use global,  only : zz


      implicit none

      external EOSINI, EOS, OPAINI_ROS, OPACITY_ROS, OPAINI_PLK, OPACITY_PLK

      real(kind=rkind) :: dOPAdTgas(np)
      logical, save    :: init = .true.


      if (init) then
         write(66,"(a)") "matfct.f90:   verwende reale Zustandsfunktionen, zentr. Matfkt durch Mittelung"
         init = .false.

         call EOSINI(EOSfile)
         call OPAINI_ROS(OPArosfile)
         if (Planck) call OPAINI_PLK(OPAplkfile)
      end if




!-----------------------------------------------------------------------
!    Aufruf der Zustandsgleichung  --  neuer Zeitpunkt
!-----------------------------------------------------------------------

      call EOS("PRE","NEW",np,rho0(1:np),X(ME,1:np),Pgas0(1:np),dPgas0dD(1:np),dPgas0dE(1:np) )
      call EOS("TEM","OLD",np,rho0(1:np),X(ME,1:np),Tgas(1:np) ,dTgasdD(1:np) ,dTgasdE(1:np)  )


!-----------------------------------------------------------------------
!    Berechnung der Zustandsgleichung  --  zeitlich zentriert
!-----------------------------------------------------------------------

      Pgas0Z(:)     = zz * Pgas0(:) + (z1-zz) * Pgas0A(:)
      dPgas0ZdD(:)  = zz * dPgas0dD(:)
      dPgas0ZdE(:)  = zz * dPgas0dE(:)

      TgasZ(:)      = zz * Tgas(:) + (z1-zz) * TgasA(:)
      dTgasZdD(:)   = zz * dTgasdD(:)
      dTgasZdE(:)   = zz * dTgasdE(:)




!-----------------------------------------------------------------------
!    Aufruf der Rosseland Opazitaeten  --  neuer Zeitpunkt
!-----------------------------------------------------------------------

      call OPACITY_ROS(np,rho0(1:np),Tgas(1:np),OPAros(1:np),dOPArosdD(1:np),dOPAdTgas(1:np))

! Umrechnung der Ableitungen

      dOPArosdD(1:np) =  dOPArosdD(1:np) + dOPAdTgas(1:np) * dTgasdD(1:np)
      dOPArosdE(1:np) =                    dOPAdTgas(1:np) * dTgasdE(1:np)


!-----------------------------------------------------------------------
!    Berechnung der zeitlich zentrierten Rosseland Opazitaeten
!-----------------------------------------------------------------------

      OPArosZ(:)    = zz * OPAros(:) + (z1-zz) * OPArosA(:)
      dOPArosZdD(:) = zz * dOPArosdD(:)
      dOPArosZdE(:) = zz * dOPArosdE(:)




      if (Planck) then
!-----------------------------------------------------------------------
!    Aufruf der Planck Opazitaeten  --  neuer Zeitpunkt
!-----------------------------------------------------------------------

      call OPACITY_PLK(np,rho0(1:np),Tgas(1:np),OPAplk(1:np),dOPAplkdD(1:np),dOPAdTgas(1:np))

! Umrechnung der Ableitungen

      dOPAplkdD(:) =  dOPAplkdD(:) + dOPAdTgas(:) * dTgasdD(:)
      dOPAplkdE(:) =                 dOPAdTgas(:) * dTgasdE(:)


!-----------------------------------------------------------------------
!    Berechnung der zeitlich zentrierten Planck Opazitaeten
!-----------------------------------------------------------------------

      OPAplkZ(:)    = zz * OPAplk(:) + (z1-zz) * OPAplkA(:)
      dOPAplkZdD(:) = zz * dOPAplkdD(:)
      dOPAplkZdE(:) = zz * dOPAplkdE(:)


      else ! otherwise we use OPAplk = OPAros

         OPAplk(:)     = OPAros(:)
         dOPAplkdD(:)  = dOPArosdD(:)
         dOPAplkdE(:)  = dOPArosdE(:)
         OPAplkZ(:)    = OPArosZ(:)
         dOPAplkZdD(:) = dOPArosZdD(:)
         dOPAplkZdE(:) = dOPArosZdE(:)

      end if


END SUBROUTINE MATFCT





SUBROUTINE UPDATE_MATFCT
!
! Verwendungszweck:
!
!     Updaten der alten Materialfkt
!
!     'UPDATE_MATFCT' wird im code 2 mal aufgerufen: in 'getready' und in 'newton'
!
!-----------------------------------------------------------------------

      implicit none


      OPArosA(:) = OPAros(:)
      OPAplkA(:) = OPAplk(:)
      TgasA(:)   = Tgas(:)
      Pgas0A(:)  = Pgas0(:)

END SUBROUTINE UPDATE_MATFCT


END MODULE MATVAR
